Flowing sand - a possible physical realization of Directed Percolation 
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A simple model for flowing sand on an inclined plane is introduced. The model is related to recent 
experiments by Douady and Daerr [Nature 399, 241 (1999)] and reproduces some of the experimen- 
tally observed features. Avalanches of intermediate size appear to be compact, placing the critical 
behavior of the model into the universality class of compact directed percolation. On very large 
scales, however, the avalanches break up into several branches leading to a crossover from compact 
to ordinary directed percolation. Thus, systems of flowing granular matter on an inclined plane 
could serve as a first physical realization of directed percolation. 



PACS numbers: 45.70.Ht, 64.60.Ht, 64.60.Ak 
I. INTRODUCTION 

Directed Percolation (DP) is perhaps the simplest 
model that exhibits a non-equilibrium phase transition 
between an "active" or "wet" phase and an inactive "dry" 
one ||l[ . In the latter phase the system is in a single "ab- 
sorbing" state; once it reaches the completely dry state, 
it will always stay there. 

Interest in DP mainly stems from universality of the 
associated critical behavior. It is believed that transi- 
tions in all models with an absorbing state belong to 
the DP universality class (unless there are some special 
underlying symmetries). DP exponents were measured 
for an extremely wide variety of models. Even though 
the exponents have not yet been calculated analytically, 
their values (especially in 1-1-1 dimensions) are known 
with very high precision 

Despite the preponderance of models in the DP univer- 
sality class, so far no physical system has been found to 
exhibit DP behavior. Indeed, as noted by Grassberger, 

"...t/iere is still no experiment where the crit- 
ical behavior of DP was seen. This is a very 
strange situation in view of the vast and suc- 
cessive theoretical efforts made to understand 
it. Designing and performing such an experi- 
ment has thus top priority in my list of open 
problems" 

The purpose of this paper is to point out that a simple 
system of sand flow on an inclined plane, that has re- 
cently been introduced and studied by Daerr and Douady 
(DD), may well be the first physical realization of a tran- 
sition in the DP universality class jJJ^. In Sec. II we 
describe these experiments in fair detail. The data pre- 
sented by DD is of qualitative value and raises serious 
questions regarding the applicability of DP. In particular, 
the observed shapes of wet clusters differ from those seen 
in standard DP simulations; they are much more com- 



pact. Since the corresponding model, called Compact 
Directed Percolation (CDP), is unstable against pertur- 
bations towards the standard DP behavior [|| , the latter 
is the generic case expected to occur (if no parameters 
were fine-tuned to place the system in the CDP class). 

This motivated us to look for a simple model which 
is defined in terms of dynamic rules that can plausibly 
be related to the experiments and, at the same time, ex- 
hibit features that look like the experimentally observed 
ones. Whether the transition exhibited by such a model 
does belong to the DP universality class remains to be 
investigated. 

Such a model is introduced in Sec. III. It is a directed 
sandpile model, which is simpler than the one introduced 
and analyzed by Tadic and Dhar Q]; here the system is 
reset to a uniform initial state after each avalanche. In 
Sec. IV we show the outcome of some simulations. The 
avalanches (observed in the active phase) reproduce the 
experimental observations quite well. We establish the 
existence of a transition from an active to an inactive 
phase. However, the critical behavior extracted from 
these figures does not seem to be in the DP universal- 
ity class, rather, it seems close to CDP. As it turns out, 
this CDP type critical behavior is only a transient: the 
true critical behavior is of the DP type, but can only be 
seen after a very long crossover regime, in which the ex- 
ponents are those of CDP. This observation is based on 
a careful numerical study, which is presented in Sees. V 
and VI. 

Our conclusion is that the DD experiment does serve 
as a possible realization of a DP-type transition. Obser- 
vation of DP exponents may be tricky as a substantial 
crossover regime may mask the true critical behavior, and 
one should try to find methods to shorten this regime. 

Finally we should note that the DD system is a simple 
case of Self Organized Criticality (SOC). Without any 
fine tuning, the system "prepares itself" at the critical 
point of a DP type transition. The way in which this 
happens differs from standard SOC models |§] in which 
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FIG. 1. Schematic stability diagram of 

the DD-experiment. A layer of thickness h is dynamically 
stable below a certain threshold hdif) (solid line). Due to 
friction forces non-moving layers remain stable in the mixed 
region below static stability limit hs{^p) (dashed line). In the 
present work we investigate the properties in the vicinity of 
the dynamic phase transition line, as indicated by the arrows. 



a slow driving force (acting on a much time scale smaller 
than that of the system's dynamic response) causes evo- 
lution to a critical state. In the present case avalanches 
are started by hand one by one. 



II. THE DOUADY-DAERR EXPERIMENT 

The experimental apparatus consists of an inclined 
plane (size of about Im) covered by a rough velvet cloth; 
the angle of inclination ipQ can be varied. Glass beads 
(e.g. "sand") of diameter 250-425 ^m |^ are poured uni- 
formly at the top of the plane and flow down while a thin 
layer of thickness h — hd{'^o), consisting of several mono- 
layers, settles and remains immobile. At this thickness 
the sand is dynamically stable] the value of hd decreases 
with an increasing angle of inclination. 

For each Lp^ there exists another thickness hg with 
hsispo) > hd{fQ), beyond which a static layer becomes 
unstable. Hence there exists a region (see Fig. |l|) in the 
((/?, h) plane, in which a static layer is stable but a flowing 
one is unstable. We can now take the system, that settled 
at hd{(po), and increase its angle of inclination to stay- 
ing within this region of mixed stability. The layer will 
not flow spontaneously, but if we disturb it at the top, 
generating a flow near the perturbation, the flow will per- 
sist and an avalanche will be generated, leaving behind a 
layer of thickness hd{(p). These avalanches had the shape 
of a fairly regular triangle, with opening angle 9. As the 
increment of the inclination 



Aip = ip - (fio 

decreases, the value of 0{A(p) decreases as well and the 
area affected by the avalanche decreases, vanishing as 
A(p 0. This calls for testing a power law behavior of 
the form 



(1) 



If instead of increasing ip we lower the plane, i.e., go to 
Aip < 0, our system, whose thickness is hd{ipo), is below 
the present thickness of dynamic stability, hd{(p). We 
believe that in this case an initial perturbation will not 
propagate, it will rather die out after a certain time (or 
beyond a certain size ^|| of the transient avalanche). As 
the deviation \A(p\ decreases, we expect the size of the 
transient active region to increase, i.e., the decay length 
should grow according to a power law 



(2) 



Hence, by pouring sand at inclination ipo, DD produced 
a self-organized critical system. The system is precisely 
at the borderline (with respect to changing the angle) 
between a stable regime < (po in which perturbations 
die out and an unstable one, ip > ipo, where perturbations 
persist and spread. 

Once this connection has been made, it is natural to 
associate this system with the problem of DP. Denote 
by p either the site or bond percolation probability and 
by pc its critical value (i.e., for p > pc the system is in 
the active phase). We associate the change in tilt with 
p — Pc, assuming that near the angle of preparation the 
behavior of the sand system is related to a DP problem 
with 



Alp = ip - ipo (xp- Pc 



(3) 



Hence, the exponent should be compared with the 
known values for DP and CDP. The exponent x in Eq. (|l|) 
can also be measured and compared with 

tane^U/^W-i^^p-"^ ■ (4) 



III. THE MODEL 

Our aim is to write down a simple model based on the 
physics of flowing sand. We adopt the observation made 
by DD, that in the regime of interest (i.e., for tilt angles 
close to po) grains of the top layer of sand rest on grains 
of the layer below (rather than on other grains of the top 
layer )|^. Hence the lower layers provide for the top one a 
kind of washboard potential, as depicted in Fig. 



^ This holds for (p < (fio and also for p> > ipo, as long as we 
stay within the region of mixed stability. 
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FIG. 2. The top layer of sand in an effective washboard 
potential. 

We further assume that only the top layer participates 
in an avalanche and therefore place the grains of this 
layer on the sites of a regular square lattic^ (see Fig. ||) . 
At any given time a particular horizontal row of grains 
may become active, while at the next time step the activ- 
ity may be transferred to the row beneath. The physical 
picture that underlies the model is as follows. A grain 
G may become active if at least one of the neighboring 
grains in the row above it has been active at the previous 
time step. These grains may then transfer energy to G; 
if AE{G), the total energy transferred to G, exceeds the 
barrier Ef, of the washboard, G becomes active. An ac- 
tive grain "rolls down" at the next time step and collides 
with the grains of the next row. The energy it brings to 
these collisions is 1 + AE{G), where 1 is the potential 
energy due to the height difference between two consec- 
utive rows. A fraction / of its total energy is dissipated, 
while the rest is divided stochastically among its three 
neighbors from the lower row. 

The model is hence defined in terms of two variables; 
an activation variable, 

gt _ \ ^ if grain {t, i) active, 
* 1 otherwise, 

and an energy variable E^. The index t denotes rows of 
our square lattice and time; at time t we update the 
states of the grains belonging to row t. Energy is mea- 
sured in units of the difference between two successive 
minima of the potential (see Fig. ||). The model is con- 
trolled by two parameters, namely 

Eb , the barrier height, and 

/ , the fraction of dissipated energy. 

The dynamic rules of our model are defined in terms 
of these variables and parameters as follows. For given 
values of activities S* and energies Ef we first calculate 
the energy transferred to the grains of the next row t + 1. 
To this end we generate for each active site Sf = 1 three 
random numbers, zl(6) (with 6 = ±1, 0) in a way that 



■^We chose to work with a square lattice, but could have used 
a triangular one as well, with each site communicating with 
two neighbors above and two below. 



Eo 




FIG. 3. Energy transfer between grains on a square lat- 
tice. 

E = (5) 

s=±i,o 

The energy transferred to grain {t-\-l,i) is then given by 
AEl+' = (1-/) E SlsEl,zl,iS). (6) 

<5=±1,0 

The values of these energies determine the activation of 
the grains of row t + 1: 

ot+i _ f 1 if AEl+' > Eb , , . 

" \ if AEl+'<Eb. 

The energies of the active grains are set according to 

El+^ = (l + AS*+i). (8) 

The meaning of these rules, in words, is obvious: the 
energy of site i at time i + 1 is obtained by identifying, 
among its three neighbors of the preceding row, those 
sites (or grains) that were active at time t. At each such 
active site (t, i) we generated three random numbers zj(S) 
which represent the fraction of energy transferred from 
the grain at site {t,i) to the one at {t + l,i + S). We 
add up the energy contributions from these active sites; 
the fraction 1 — / is not dissipated and compared to the 
barrier height Eb- If the acquired energy AE^'^^ exceeds 
Eb , site {t + becomes active, rolls over the barrier 
bringing to the collisions (at time t + 2) the acquired 
energy calculated above and its excess potential energy 
(of value 1). 



IV. SHORT-TIME SIMULATIONS AND 
QUALITATIVE DISCUSSION OF THE 
TRANSITION 

Let us consider the behavior of our model as we vary 
Eb at a fixed value of the dissipation. We expect that 



3 



\ \ 
\ \ 
\ \ 

\ \ 




\ 
\ 

V ^ 


stable 


mixed 


\ 

\ 



0.2 0.4 0.6 0.8 1 

f 

FIG. 4. Phase diagram of the model for flowing sand. The 
full line represents the phase transition line. The mean-field 
approximations of Eqs. (^o|) and ( |l3| ) arc shown as dotted and 
dashed lines, respectively. 



for small values of Ei, an active grain will activate the 
grains below with high probability; avalanches will prop- 
agate downhill and also spread sideways. For a strongly 
localized initial activation we should, therefore, observe 
activated regions of triangular shape. As Ej, increases, 
the rate of activation decreases and the opening angle 9 
of these triangles should decrease, until Ei, reaches a crit- 
ical value E^, beyond which initial activations die out in 
a finite number of time steps (or rows). These expecta- 
tions are indeed borne out by simulations of the model: 
the critical value E^ depends on the dissipation / and 
the resulting phase transition line is shown in Fig. as a 
solid line. 

In order to understand this transition qualitatively, let 
us consider a simple mean-field type approximation, in 
which all stochastic variables are replaced by their aver- 
age values. 

We consider an edge separating an active region from 
an inactive one at time t: sites to the left of i and i itself 
are wet, whereas i + l,i + 2, ... are dry. Will the rightmost 
wet site be wet or dry at the next time step? Assuming 
that all wet sites at time t have the same energy i?*, in 
our mean-field type estimate the energy delivered to site 
i at time i -I- 1 is 



t+i 



'-(i-m + AE*)., 



(9) 



where we set in Eq. ^ all z{S) — 1/3. At the critical 
point we expect all energies just to be sufficient to go over 
the barrier; hence set AK*+^ AE^ = Ej^ in Eq. (|). 
Solving the resulting equation yields 

E^-^^. (10) 
''1 + 2/ ^ ' 

In Fig. H this rough estimate of the transition line is 



FIG. 5. Schematic drawing of the energy profile of a com- 
pact cluster in the improved mean field approximation. 

shown as a dotted line. 

This simple calculation captures the physics of the 
problem. However, it is easy to improve it in the fol- 
lowing way. As before, we assume the energy of toppling 
grains to be distributed equally among the three neigh- 
bors of the subsequent row. However, we no longer as- 
sume all active sites to carry the same energy, instead we 
compute the energy profile at the edge of a cluster. To 
this end let us consider a semi-infinite cluster with Si ~ 1 
for i < and S"* = for i > 0. According to Eq. (||), 
we are looking for a stationary solution of the equation 
of motion 




AEl^^ if i < 
if i = 
if i > 



where AEq — E^. The corresponding stationary solution 
reads 

AEf'^' = Ebuik - Egap exp {ai) , {i < 0) (11) 
where 



Ebuik = {i-.n/.f, 

_2 + f- V12/-3/2 



E, 



gap 



a = arccosh 



2/(1-/) 
2 + / 



(12) 



2-2/ 

Thus, the critical threshold is given by the expression 



E^. 



2/' - 5/ + v/12/ ^ 3/^ 
2/(/ - 1) 



(13) 



which slightly improves the mean field result (|T^), espe- 
cially for small values of / (see dashed line in Fig. ^).The 
energy profile decreases at the edges of the cluster and 
saturates in the bulk at Ehuik, as shown in Fig. ^ 

The connection of our model to the experimental con- 
ditions is based on the assumption that the tilt angle of 
the experiment tunes the ratio between the barrier height 
and the difference of potential energies between two rows. 
If the system has been prepared at some tpo , we raise the 
tilt angle to (p; perturbing the system in this region of 
mixed stability will generate an avalanche. 

That is, for (p > (po we have Eb < E^. As the tilt 
angle is reduced, the size of Eh (measured in units of the 
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EtFO.300 Ei=0.377 

FIG. 6. Typical avalanches starting from a single seed 
with dissipation / — 0.5 far away and close to criticality. 

potential difference) increases, until it reaches its critical 
value precisely at cpQ. Thus increasing Eb in the model 
corresponds to lowering the tilt angle towards the value 
at which the system has been prepared and, as such, is 
precisely the boundary of dynamic stability. 

Hence to reproduce the experiment we were looking for 

1. fairly compact triangular regions of activation for 

2. a varying opening angle of these triangles which 
should go to zero as Eb approaches E^ from below. 

The number of "time steps" that correspond to the DD 
experiment can be estimated as the number of rows of 
beads from top to bottom of the plate, i.e. about 3000. 

We simulated the model defined in Eqs. (^)-(||) to check 
whether it is possible to reproduce the qualitative fea- 
tures of the experiment. Indeed we found this to be the 
case, as can be seen in Fig. ||. The two avalanches were 
produced for dissipation / = 0.5, activating a single site 
a,t t — 0, to which an initial energy of Eq — 500 was as- 
signed^. The avalanches were compact, triangular, and 
with fairly straight edges. The edges became rough only 
when Eb was very close to its critical value, as can be 
seen on the right hand side of Fig. ^. The opening angle 
of the active regions 9 decreased as Eb increased towards 
E^, which is shown in Fig. 0. From these simulations we 
obtain the estimate (see Eq. (^) 

x = i'\\-iy±= 0.98(5) ~ 1 . (14) 
We predict that measuring the dependence of the 



FIG. 7. Opening angle 9 of activated triangular re- 
gions as a function of the distance from criticality in a dou- 
ble-logarithmic representation. 

avalanche opening angle on Aip in the experiment should 
also give a linear law. 

Furthermore, the density of active sites in the interior 
of the triangular regions is found to be almost constant, 
indicating a first-order transition. These results suggest 
that the transition belongs to the CDF universality class, 
which is characterized by the critical exponents ^ 

i^\\=2, = /3 = 0. (15) 

These observations pose, however, a puzzle; since we be- 
lieve that DP is the generic situation, we would expect to 
find non-compact active regions and DP exponents. In 
the following Section we present a careful numerical anal- 
ysis of the critical behavior of our model which resolves 
this problem: the exponents seen in our simulations (and 
in the experiment) should cross over to the DP values, 
but only if one gets very deep into the critical region. 



V. CROSSOVER TO DIRECTED PERCOLATION 

The linear law observed in Fig. can be explained by 
assuming compact clusters whose temporal evolution is 
determined by the fluctuations of their boundaries. The 
boundaries perform an effective random walk with a spa- 
tial bias proportional to Eb — E^. Therefore, the critical 
model should behave in the same way as a Glauber-Ising 
model at zero temperature, i.e., the transition should be- 
long to the GDP universality class. However, according 
to the DP conjecture [|l^ any continuous spreading tran- 
sition from a fluctuating active phase into a single frozen 
state should belong to the universality class of directed 



has been dissipated. 

Note that after less than 20 time steps all the initial energy 
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FIG. 8. Mean survival probability P{t) of the toppling 
process at criticality averaged over 50 000 independent runs. 
The predicted slopes for CDP and DP are indicated by dotted 
and dashed lines, respectively. 
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FIG. 9. Average number of active sites N{t). The ex- 
pected slopes are indicated in the same way as in Fig. pi 



percolation (DP), provided that the model is defined by 
short range interactions without exceptional properties 
such as higher symmetries or quenched randomness (see 
Sec. VI). Clearly, the present model fulfills these require- 
ments. It has indeed a fluctuating active state and ex- 
hibits a phase transition into a single absorbing state 
which is characterized by a positive one-component or- 
der parameter. According to these arguments, the phase 
transition should belong to the DP universality class. 

In order to understand this apparent paradox we per- 
form high-precision Monte-Carlo simulations for dissi- 
pation/ = 0.5. We employ time-dependent simula- 
tions ||ll|], i.e., we topple a single grain in the center and 
analyze the properties of the resulting cluster. As usual 
for this type of simulations, we measure the survival prob- 
ability P{t), the number of active sites N{t), and the 
mean square spreading from the origin B?{t) averaged 
over the surviving runs. At criticality, these quantities 
are expected to show an asymptotic power law behavior 
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FIG. 10. Mean square spreading from the origin averaged 
over surviving runs. In order to demonstrate the crossover 
from CDP to DP we divided R^{t) by t. The expected slopes 
(2 — z)/z are indicated by dotted and dashed lines. 



P(t) ~ t 



N{t) - 



i?2(t) ^ t 



II 7. 



(16) 



where 5, and z are critical exponents which label the 
universality class. In the case of CDP these exponents 
are given by 

z = 2, (17) 



J = l/2, r; = 0. 



whereas DP is characterized by the exponents 

(5 = 0.1595, 77 = 0.3137, z= 1.5807. (18) 

In order to eliminate finite-size effects, we use a dy- 
namically generated lattice adjusted to the actual size 
of the cluster. Moreover, we observe that the initial 
non-universal transient is minimal if an excitation energy 
i?o — 15 is used. Detecting deviations from power-law 
behavior in the long-time limit we estimate the critical 
energy by — 0.385997(5). Our numerical results (ob- 
tained from simulations at the critical point) are shown 
in Figs. |^|0|. In all measurements we observe different 
temporal regimes: 

1. During the first few time steps, the activation en- 
ergy is distributed to the nearest neighbors whereby 
the cluster grows at maximal speed. Therefore, the 
survival probability P(i) is 1 and the particle num- 
ber N(t) grows linearly. 

2. In the intermediate regime, which extends up to a 
few hundred time steps, the inactive islands within 
the cluster are not yet able to break up the cluster 
into separate parts. Thus, the cluster can be con- 
sidered as being compact and the temporal evolu- 
tion is governed by a random walk of its bound- 
aries. In this regime we observe a power-law be- 
havior with CDP exponents (indicated by dotted 
fines in Figs. HO). 



6 



CDP: 



o 
o 
o 





100 sites 800 sites 

FIG. 11. Typical clusters generated at criticality on small 
and large scales, illustrating the crossover from CDP to DP. 



3. The intermediate regime is followed by a long 
crossover from CDP to DP extending over almost 
two decades up to more than 10'* time steps^ 

4. Finally the system enters an asymptotic DP regime 
(indicated by dashed lines in Figs. 



The crossover from CDP to DP is illustrated in Fig. |T^. 
Two avalanches are plotted on different scales. The left 
one represents a typical avalanche within the first few 
thousand time steps. As can be seen, the cluster appears 
to be compact on a lateral scale up to 100 lattice sites. 
However, as shown in the right panel of Fig. |ll], after a 
very long time the cluster breaks up into several branches. 
The right hand figure shows a typical cluster on a scale of 
150 000 time steps, where the branches still have a certain 
characteristic thickness. Going to even larger scales the 
width of the branches becomes irrelevant and we obtain 
the typical patterns of critical DP clusters. 

In comparison with ordinary DP lattice models, in the 
present model the observed crossover is unusually slow. 
This due to short-range correlations between active sites 
leading to active branches with a certain typical thickness 
^act- In ordinary DP lattice models the average size of 
active branches is of the order of a few lattice spacings. 
In the present case, however, we find a much larger value 
Uct « 20. 

Based on this observation, the typical crossover time tc 
can be approximated as follows. In order to cross over to 
DP, the average size of inactive regions between neighbor- 
ing branches ^inact has to become larger than the thick- 
ness of the branches ^act- In Fig. |l^ we plot both quan- 
tities as a function of time at criticality, using a lattice 
with N — 2^^ sites and homogeneous initial conditions 
£1=° = 2. Initially ^act = N and ^inact = 0. As time 




t 

FIG. 12. Mean sizes of active and inactive regions as a 
function of time, starting from homogeneous initial conditions 
with dissipation / = 0.5 (see text). The inset shows the 
saturation value of (^act as a function of the dissipation /. 



evolves, the average size of active branches decreases and 
saturates at a constant value ^act ~ 20. However, the 
average size of inactive regions £,inact continues to grow 
and exceeds ^act at time tc « 10^. As can be seen, this 
provides a good estimate of the typical time where the 
critical behavior of the system crosses over to DP. 

In order to observe the crossover experimentally, it 
would be interesting to know how the crossover time tc 
can be reduced. To this end we measure £,act for sev- 
eral values of the dissipation / (see inset of Fig. ^l]). It 
turns out that by increasing / the typical size of active 
branches can be decreased down to 10 lattice spacings. 
Consequently, the crossover time can be reduced by more 
than one decade. Hence, for an experimental verification 
of DP, systems with high dissipation are more appropri- 
ate. 

The influence of the dissipation can easily be ex- 
plained within the improved mean field approximation 
of Sect. IV. Clearly, the stability of a cluster against 
breakup into several branches by fluctuations depends on 



the energy gap Eg 



El 



bulk 



Ec- As can easily be ver- 



ified, this energy difference (and therewith the stability 
of compact clusters) decreases with increasing dissipation 
/, explaining the observed /-dependence. 



VI. THE EFFECT OF RANDOMNESS 



*Note that the crossover in the present model is different 
from the one studied in |l^, where inhomogeneous interac- 
tions at the cluster's boundaries were assumed. 



The above model describes the physics of flowing sand 
in a highly idealized manner. In particular, it ignores the 
fact that spreading avalanches may be subjected to frozen 
disorder. For example, irregularities of the plate and the 
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velvet cloth could lead to quenched randomness in the 
equations of motion. Moreover, the system prepares itself 
in an initial state which is not fully homogeneous. Thus, 
we have to address the question to what extent quenched 
randomness will affect the expected crossover to DP. 

Certain types of quenched disorder are known to 
change the critical behavior of DP. For example, Mor- 
eira and Dickman studied the diluted contact process 
with spatially quenched disorder Even for small 

amplitudes quenched randomness was found to destroy 
the DP transition, turning algebraic into logarithmic 
laws. Janssen |14 confirmed and substantiated these 
findings by a field-theoretic analysis. Recently Cafiero 
et al. ||l^ mapped DP with spatially quenched disorder 
onto a non-markovian process with memory exhibiting 
the same nonuniversal properties. The memory is due 
to the formation of bound states of particles in those re- 
gions where the percolation probability is very high. As 
shown by Webman et al., these bound states give rise to 
a glassy phase separating active and inactive parts of the 
phase diagram Pq| . Similar nonuniversal properties were 
also observed in DP processes with temporally quenched 
disorder iQ. 

In all cases investigated so far, quenched disorder de- 
stroys the DP transition. However, the disorder in the 
DD experiment is different in nature. Clearly, it is nei- 
ther spatially nor temporally quenched, rather it depends 
on both space and time. On the level of our model we 
may think of randomly varying energy barriers 



Eb ^ Eb + Aij{x,t) , 



(19) 



where the amplitude A controls the intensity of disorder. 
Here r]{x,t) is a white Gaussian noise specified by the 
correlations 



rj{x, t)rj{x',t') = S'^ix - x')6{t - t') , 



(20) 



where d = 1 denotes the spatial dimension. In the stan- 
dard situation of quenched noise of this type r]{x,t) is 
kept fixed while the experiment is repeated and the quan- 
tities under investigation are averaged over many inde- 
pendent avalanches. Yet in the DD experiment, the sit- 
uation is different. Here once the sand has been poured, 
a particular realization of the random variables has been 
selected. However, there is no process to repeat the ex- 
periment over and over again with a fixed ry(x, t). Rather, 
after each avalanche the system is prepared again (by 
pouring sand or by starting an avalanche elsewhere). 
Hence the averaging process is done simultaneously over 
the 77(2;, t) and the stochastic dynamic process that gen- 
erates the avalanches. This type of averaging is of the 
annealed type and therefore less likely to alter the criti- 
cal behavior than its quenched version. 

In order to find out whether fully quenched disorder 
affects the asymptotic critical behavior of DP, we simu- 
lated a directed bond percolation process with randomly 
distributed bond probabilities between p* and 1. For 



p* = 0.289(1), we find asymptotic power laws with DP 
exponents, indicating that the transition is not affected 
by spatio-temporally quenched noise. Therefore, we ex- 
pect the same to be true in the case of annealed disorder 
in our model for flowing sand. 

To support this point of view, we study the case of 
quenched randomness in the DP Langevin equation pq| 



dtp{x, t) = ap{x, t) - gp^{x, t) + 

DVp(x, t) + V^p{x,t)i{x, t) , 



(21) 



where p{x, t) is the particle density and a represents the 
percolation probability. £,{x,t) is a Gaussian white noise 
which represents the intrinsic randomness of the DP pro- 
cess. At the critical dimension d = 4, where fluctuations 



start to contribute, the Langevin equation (21) is invari- 
ant under scaling transformations x — > bx, t — > iP't, and 
p^h-'^p. 

In order to include spatio-temporally quenched ran- 
domness, we allow for small variations of a, i.e., we add 
the term 

Ap{x,t)T]{x,t) 

on the right hand side of Eq. (pl|). However, as can be 
shown by simple dimensional analysis, this term is irrele- 
vant m. d — A dimensions, i.e., it decreases and eventually 
vanishes under scaling transformations. This observation 
strongly supports the result that the DP transition in our 
model is indeed not affected by quenched randomness. 

We emphasize that the irrelevance of quenched ran- 
domness in our model is due to the special role of 'time' 
which coincides with the vertical coordinate of the plane. 
That is, for each time step the stochastic processes take 
place in a different random environment. To that ex- 
tent the DD experiment differs from other DP-related 
experiments such as catalytic reactions where spatially 
quenched disorder affects the critical behavior. 



VII. CONCLUSIONS 

We introduced a simple model for flowing sand on an 
inclined plane. The model is inspired by recent exper- 
iments and reproduces some of the observed features. 
In contrast to the experiment, which prepares itself in 
a self-organized critical state, our model needs to be 
tuned to a critical point by varying the energy barrier 
Eb- At criticality the system undergoes a nonequilibrium 
phase transition from an inactive (dry) phase with finite 
avalanches to an active (wet) phase where the mean size 
of avalanches diverges. Analyzing the critical behavior 
near the transition, we obtained the following results: 

1. On short scales, i.e., on scales considered in the DD 
experiment, the model reproduces the experimen- 
tally observed triangular compact avalanches. In 
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the active phase their opening angle 9 is predicted 
to vary hnearly^ with Aip. 

2. On very large scales the critical behavior of the 
model crosses over to ordinary DP. Thus, the DD 
experiment could serve as a first physical real- 
ization of directed percolation. Crossover to DP 
is seen in the model after about 10'' time steps, 
whereas the DD experiment stops at about 3000 
steps (i.e. rows of beads). Hence in order to ob- 
serve the crossover in the experiment, larger system 
sizes and/or smaller beads would be required. 

3. We have shown that quenched randomness with 
short-range correlations due to irregularities in the 
experiment should not affect the asymptotic critical 
behavior. 

4. The typical time needed to cross over to DP is 
found to decrease with increasing dissipation. 

Thus, in order to create experimental conditions fa- 
voring a crossover to DP, we suggest to use small glass 
beads, large system sizes, and an initial angle ipQ where 
the dissipation of energy per toppling grain is maximal. 
For physical reasons we would expect the dissipation to 
be maximal for small angles (^o i but this has to be verified 
in the actual experiment. 

As a necessary precondition for a crossover to DP, 
compact clusters must be able to split up into several 
branches, as illustrated in Fig. 0. Thus, before measur- 
ing critical exponents, this feature has to be tested ex- 
perimentally. To this end the DD experiment should be 
performed repeatedly at the critical tilt cp = ipQ. In most 
cases the avalanches will be small and compact. However, 
large avalanches, reaching the bottom of the plate, will 
sometimes be generated. If these avalanches are non- 
compact (consisting of several branches) we expect the 
asymptotic critical behavior to be described by DP. Only 
then it is worthwhile to optimize the experimental setup 
and to measure the critical exponents quantitatively. 
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^Note added after submission: This prediction has to be 
compared with the model proposed by Bouchaud et al, ^] 
which predicts the exponent x = 1/2. 
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